# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")

pacman::p_load(tidyverse, # data wrangling + plotting
               gt, # to create tables
               brms, # to fit Bayesian models
               here, # reproducible file paths
               rio, # to import files
               ggdist, # to plot distributions 
               tidybayes, # to plot distributions 
               PNWColors,
               flextable, # table
               patchwork) # arrange plots

# Load data
d_logOR = readRDS(
  here("03_updated_analyses_who", "output", "data", "effect_sizes.Rds")
  )

# Load main model
m1 = readRDS(here("03_updated_analyses_who", "output", "fits", "m1.Rds"))

# Load function
source(here("03_updated_analyses_who/functions/diag_plot.R"))

In this file, we will discuss the rationale and fit models in which we use simulated randomized controlled trials as priors.

Details

The main Bayesian meta-regression random-effect model is:

\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu & = \alpha + \beta_{SOO} x_i + \beta_{NIV}x_i\\ \\ \alpha & \sim \operatorname{Normal}(0, 0.82^2) \tag{Priors} \\ \beta_{SOO}, \beta_{NIV} & \sim \operatorname{Normal}(0, 1.5^2) \\ \tau & \sim \operatorname{Half-Normal}(0.5^2) \\ \end{align*} \]

where \(y_i\) is the observed mean log odds ratio of tocilizumab versus control and \(\sigma_i^2\) is the known sampling variance in study \(i\). \(\alpha\) is the intercept, which represents the overall effect of tocilizumab in patients on invasive mechanical ventilation. \(\beta_{SOO}\) is the difference between patients on simple oxygen only (SOO) and patients on invasive mechanical ventilation (\(\alpha\)). \(\beta_{NIV}\) is the difference between patients on noninvasive ventilation (NIV) and patients on invasive mechanical ventilation (\(\alpha\)). Both coefficients are multiplied by \(x\), which is dummy-coded. Lastly, \(\tau^2\) represents the between-study heterogeneity.

### Overall mean odds ratio in WHO's meta-analysis on tocilizumab in patients
# using corticosteroids:

# The mean odds ratio in patients using corticosteroids (toci vs control)
# was 0.77 (page 14 in Supplement 2 doi:10.1001/jama.2021.11330)

mean_OR_WHO = 0.77 

Here, we will fit 6 models using different priors for \(\alpha\). Each prior will be tailored to simulate a randomized controlled trial (RCT) of tocilizumab versus standard of care (control). We want to simulate RCTs with the following total sample sizes: \(200, 500, 1000, 1500, 2000, 4000\).

All priors will be normally distributed and will have a mean of -0.26. This value is the equal to log(0.77), which was chosen based on WHO’s meta-analysis (page 14 in Supplement 2). This is the mean odds ratio of tocilizumab vs. control in patients using corticosteroids (overall results).

To calculate the standard deviation of each prior, we need to define certain aspects of these “simulated randomized controlled trials”:

control_risk = 0.43

# Now, we need to calculate the risk in the tocilizumab arm based on the odds ratio
# and control risk mentioned above

# Reference for formula: Box 1 in https://doi.org/10.1016/j.jclinepi.2020.08.019

toci_risk = (control_risk*mean_OR_WHO)/(control_risk*(mean_OR_WHO - 1) + 1)
  1. We assume equal allocation in both treatment arms.

  2. Adapting from the suggestions in the GRADE guidelines, we found a striking discrepancy between the control mortality risk in this data (52%) in comparison to another previously published meta-analysis (34% in patients on IVM and using corticosteroids). Thus, we have decided to use 43% (arithmetic mean between 34 and 52) as our reference risk in the IVM subgroup..

  3. The mortality risk in the tocilizumab was calculated using the following formula (Box 1 in Doi et al., 2020):

\[R_{T} = \frac{R_{C} O R}{R_{C}\left(O R-1\right)+1}\]

where \(R_{T}\) is the mortality risk in the tocilizumab group, \(R_{C}\) is the mortality risk in the control group and \(OR\) is the odds ratio mentioned above. Thus, the tocilizumab risk is equal to 0.37.

Thus, we are simulating RCTs with mean OR equal to 0.77, control risk mortality of 0.43, and tocilizumab risk of 0.37.

Based on these values, we can estimate the standard deviation (\(SD\)) with the following formula:

\[SD=\sqrt{\frac{1}{a+\frac{1}{2}}+\frac{1}{b+\frac{1}{2}}+\frac{1}{c+\frac{1}{2}}+\frac{1}{d+\frac{1}{2}}}\]

where \(a\), \(b\), \(c\) and \(d\) are number of events and follow this 2x2 table:

tribble(
  
  ~"Event", ~"Tocilizumab", ~"Control",
  
  "Death", "a", "c",
  "No death", "b", "d"
) %>% 
  flextable() %>% 
  autofit()

As similarly shown in the supplementary material of Higgins and Spiegelhalter, 2002, we can estimate these values as:

  • \(a = R_{T}SS_{T}\)
  • \(b = SS_{T} - R_{T}SS_{T}\)
  • \(c = R_{C}SS_{C}\)
  • \(d = SS_{C} - R_{C}SS_{C}\)

where \(SS_{T}\) and \(SS_{C}\) are the sample sizes in the tocilizumab and control arms, respectively. As mentioned above, we assume equal allocation in both treatment arms, thus \(SS_{T} = SS_{C}\).

# Function to calculate the standard deviation of odds ratio from simulated RCT ,
# based on risk and total sample size in each treatment arm
# Here we assume both arms have the same number of total patients (total_arm)

# Adapted from:
# Appendix in https://doi.org/10.1093/ije/31.1.96
# +
# Appendix in https://doi.org/10.1016/j.jclinepi.2008.07.006

SD_fun = function(risk_toci, # Risk in tocilizumab arm
                  risk_control, # Risk in control arm
                  total_arm){ # Total sample size in each arm
  
# Tocilizumab
a = risk_toci*total_arm                     # Deaths
b = total_arm - risk_toci*total_arm         # No deaths
# Control
c = risk_control*total_arm                  # Deaths
d = total_arm - risk_control*total_arm      # No deaths

# Standard deviation =
sqrt(1/(a + 1/2) + 1/(b + 1/2) + 1/(c + 1/2) + 1/(d + 1/2)) 

                  }
# Lastly, we should apply the function above to calculate the standard deviation
# of multiple simulated RCTs will the same risk_toci and risk_control, but 
# different amounts of patients

df = 
  tibble(risk_toci = toci_risk,
         risk_control = control_risk,
         total_arm = c(100, 250, 500, 750, 1000, 2000)) %>% 
  mutate(SD = SD_fun(risk_toci, risk_control, total_arm))

Finally, we can estimate the \(SD\) based on the 6 different sample sizes mentioned above:

df %>% 
  summarise("Sample size in each treatment arm" = total_arm,
            "Total sample size" = 2*total_arm,
            " SD" = round(SD,2)) %>% 
  flextable() %>% 
  autofit()

In summary, there are 4 parameters in these models, which are \(\alpha\), \(\beta_{SOO}\), \(\beta_{NIV}\), and \(\tau\). The priors for the latter three parameters will be the same in every model:

\[ \begin{align*} \beta_{SOO}, \beta_{NIV} & \sim \operatorname{Normal}(0, 1.5^2)\\ \tau & \sim \operatorname{Half-Normal}(0.5^2) \end{align*} \]

On the other hand, the priors for \(\alpha\) are described as:

\[ \begin{align*} \alpha & \sim \operatorname{Normal}(\mu, SD^2) \end{align*} \]

where \(\mu\) is equal to -0.26 in every model, but \(SD\) range from 0.29 to 0.06, as described in the table above.

Let’s fit these 6 new models:

# Priors

intercept100 = prior_string(
  paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[1,4]], # SD, total_arm = 100
  ")"),
  class = "b", coef = "Intercept"
  )

# b_Intercept ~ normal(-0.261364764134408,0.28303369902181)

intercept250 = prior_string(
  paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[2,4]], # SD, total_arm = 250
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.179552672191073)


intercept500 = prior_string(paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[3,4]], # SD, total_arm = 500
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.127092510741972)

intercept750 = prior_string(paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[4,4]], # SD, total_arm = 750
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.103805945961058)

intercept1000 = prior_string(paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[5,4]], # SD, total_arm = 1000
  ")"),
  class = "b", coef = "Intercept")

#b_Intercept ~ normal(-0.261364764134408,0.0899139032236548)

intercept2000 = prior_string(paste0(
  "normal(",
  log(mean_OR_WHO), # Mean, log(0.77)
  ",",
  df[[6,4]], # SD, total_arm = 2000
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.0635949873256457)

SOO = prior(normal(0, 1.5), class = "b", coef = "oxygenlow") # beta_Simple oxygen only

NIV = prior(normal(0, 1.5), class = "b", coef = "oxygenNIV") # beta_NIV

tau_study = prior(normal(0, 0.5), class = "sd") # tau_study

priors = list(
  
  c(intercept100, SOO, NIV, tau_study),
  c(intercept250, SOO, NIV, tau_study),
  c(intercept500, SOO, NIV, tau_study),
  c(intercept750, SOO, NIV, tau_study),
  c(intercept1000, SOO, NIV, tau_study),
  c(intercept2000, SOO, NIV, tau_study)
)

# Create list to store fits
# new_models = list()
# 
# # Store the main model
# 
# new_models[[1]] = m1
# 
## Supressed because we saved the models already
#
## Run the loop, to fit other models with different 
#
# for (i in 1:nrow(df)) {
# 
#   # Show what model is running
#   print(i)
#   
#   # Skip first index because main model (m1) is already stored there
#   k = i + 1
#   
#   # Re-fit the main model, but now changing the prior
#   new_fit = update(m1, 
#                prior = priors[[i]], 
#                
#                cores = parallel::detectCores(),
#                sample_prior = T,
#                chains = 4,
#                control = list(adapt_delta = .99),
#                backend = "cmdstanr",
#                seed = 123)
#   
#   # Save models
#   new_models[[k]] = new_fit
# }
# # 
# saveRDS(new_models,
#          here("03_updated_analyses_who", "output", "fits", "simulated_priors_models.rds"))

new_models = readRDS(here("03_updated_analyses_who",
                          "output",
                          "fits",
                          "simulated_priors_models.rds"))

Diagnostic plots

# Run loop
for (i in 1:(nrow(df))) {
  
  
# Skip first index (main model)
k = i + 1

  # Run custom function diag_plot, 1 per model
p = diag_plot(model = new_models[[k]],
          pars_list = c("b_Intercept", "b_oxygenlow", "b_oxygenNIV", "sd_study__Intercept"),
          ncol_trace = 4)
# Display plot
print(p)  

# Add legend to show respective priors
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), 0.25, .03, gp=grid::gpar(cex=1.3))
}

Posterior predictive check

# Run lopp
for (i in 1:(nrow(df))) {
  
# Skip first index (main model)
k = i + 1

p = pp_check(new_models[[k]], ndraws = 20)

print(p)

lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), .7, .9, gp=grid::gpar(cex=1.2))
}

Sensitivity analysis model

We will fit models assuming a priors for \(\alpha\) with the mean equal to \(log(1)\).

# Priors

intercept100 = prior_string(
  paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[1,4]], # SD, total_arm = 100
  ")"),
  class = "b", coef = "Intercept"
  )

# b_Intercept ~ normal(0,0.28303369902181)

intercept250 = prior_string(
  paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[2,4]], # SD, total_arm = 250
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.179552672191073)


intercept500 = prior_string(paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[3,4]], # SD, total_arm = 500
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.127092510741972)

intercept750 = prior_string(paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[4,4]], # SD, total_arm = 750
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.103805945961058)

intercept1000 = prior_string(paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[5,4]], # SD, total_arm = 1000
  ")"),
  class = "b", coef = "Intercept")

#b_Intercept ~ normal(-0.261364764134408,0.0899139032236548)

intercept2000 = prior_string(paste0(
  "normal(",
  log(1), # Mean, 0
  ",",
  df[[6,4]], # SD, total_arm = 2000
  ")"),
  class = "b", coef = "Intercept")

# b_Intercept ~ normal(-0.261364764134408,0.0635949873256457)

SOO = prior(normal(0, 1.5), class = "b", coef = "oxygenlow") # beta_Simple oxygen only

NIV = prior(normal(0, 1.5), class = "b", coef = "oxygenNIV") # beta_NIV

tau_study = prior(normal(0, 0.5), class = "sd") # tau_study

priors = list(
  
  c(intercept100, SOO, NIV, tau_study),
  c(intercept250, SOO, NIV, tau_study),
  c(intercept500, SOO, NIV, tau_study),
  c(intercept750, SOO, NIV, tau_study),
  c(intercept1000, SOO, NIV, tau_study),
  c(intercept2000, SOO, NIV, tau_study)
)

# Create list to store fits
# new_models = list()
# # 
# # # Store the main model
# # 
# new_models[[1]] = m1
# # 
# # ## Supressed because we saved the models already
# # 
# # # Run the loop, to fit other models with different 
# # 
# for (i in 1:nrow(df)) {
# 
#   # Show what model is running
#   print(i)
#   
#   # Skip first index because main model (m1) is already stored there
#   k = i + 1
#   
#   # Re-fit the main model, but now changing the prior
#   new_fit = update(m1, 
#                prior = priors[[i]], 
#                
#                cores = parallel::detectCores(),
#                sample_prior = T,
#                chains = 4,
#                control = list(adapt_delta = .99),
#                backend = "cmdstanr",
#                seed = 123)
#   
#   # Save models
#   new_models[[k]] = new_fit
# }
# # 
# saveRDS(new_models,
#          here("03_updated_analyses_who", "output", "fits",
#               "sensitivity_simulated_priors_models.rds"))

new_models = readRDS(here("03_updated_analyses_who",
                          "output",
                          "fits",
                          "sensitivity_simulated_priors_models.rds"))

Diagnostic plots

# Run loop
for (i in 1:(nrow(df))) {
  
# Skip first index (main model)
k = i + 1

  # Run custom function diag_plot, 1 per model
p = diag_plot(model = new_models[[k]],
          pars_list = c("b_Intercept", "b_oxygenlow", "b_oxygenNIV", "sd_study__Intercept"),
          ncol_trace = 4)
# Display plot
print(p)  

# Add legend to show respective priors
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), 0.25, .03, gp=grid::gpar(cex=1.3))
}

Posterior predictive check

# Run lopp
for (i in 1:(nrow(df))) {
  
# Skip first index (main model)
k = i + 1

p = pp_check(new_models[[k]], ndraws = 20)

print(p)

lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), .7, .9, gp=grid::gpar(cex=1.2))
}

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
##  [1] bayesplot_1.8.1 cowplot_1.1.1   patchwork_1.1.1 flextable_0.6.8
##  [5] PNWColors_0.1.0 tidybayes_3.0.1 ggdist_3.0.0    rio_0.5.27     
##  [9] here_1.0.1      brms_2.16.1     Rcpp_1.0.7      gt_0.3.1       
## [13] forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4    
## [17] readr_2.0.1     tidyr_1.1.3     tibble_3.1.4    ggplot2_3.3.5  
## [21] tidyverse_1.3.1 pacman_0.5.1   
## 
## loaded via a namespace (and not attached):
##   [1] uuid_0.1-4           readxl_1.3.1         backports_1.2.1     
##   [4] systemfonts_1.0.2    plyr_1.8.6           igraph_1.2.6        
##   [7] svUnit_1.0.6         splines_4.0.4        crosstalk_1.1.1     
##  [10] rstantools_2.1.1     inline_0.3.19        digest_0.6.27       
##  [13] htmltools_0.5.2      rsconnect_0.8.24     fansi_0.5.0         
##  [16] magrittr_2.0.1       checkmate_2.0.0      openxlsx_4.2.4      
##  [19] tzdb_0.1.2           modelr_0.1.8         RcppParallel_5.1.4  
##  [22] matrixStats_0.60.1   officer_0.4.0        xts_0.12.1          
##  [25] prettyunits_1.1.1    colorspace_2.0-2     rvest_1.0.1         
##  [28] haven_2.4.3          xfun_0.25            callr_3.7.0         
##  [31] crayon_1.4.1         jsonlite_1.7.2       lme4_1.1-27.1       
##  [34] zoo_1.8-9            glue_1.4.2           gtable_0.3.0        
##  [37] V8_3.4.2             distributional_0.2.2 pkgbuild_1.2.0      
##  [40] rstan_2.21.2         abind_1.4-5          scales_1.1.1        
##  [43] mvtnorm_1.1-2        DBI_1.1.1            miniUI_0.1.1.1      
##  [46] xtable_1.8-4         foreign_0.8-81       stats4_4.0.4        
##  [49] StanHeaders_2.21.0-7 DT_0.19              htmlwidgets_1.5.4   
##  [52] httr_1.4.2           threejs_0.3.3        arrayhelpers_1.1-0  
##  [55] posterior_1.1.0      ellipsis_0.3.2       pkgconfig_2.0.3     
##  [58] loo_2.4.1            farver_2.1.0         sass_0.4.0          
##  [61] dbplyr_2.1.1         utf8_1.2.2           labeling_0.4.2      
##  [64] tidyselect_1.1.1     rlang_0.4.11         reshape2_1.4.4      
##  [67] later_1.3.0          munsell_0.5.0        cellranger_1.1.0    
##  [70] tools_4.0.4          cli_3.0.1            generics_0.1.0      
##  [73] broom_0.7.9          ggridges_0.5.3       evaluate_0.14       
##  [76] fastmap_1.1.0        yaml_2.2.1           processx_3.5.2      
##  [79] knitr_1.34           fs_1.5.0             zip_2.2.0           
##  [82] nlme_3.1-152         mime_0.11            projpred_2.0.2      
##  [85] xml2_1.3.2           compiler_4.0.4       shinythemes_1.2.0   
##  [88] rstudioapi_0.13      gamm4_0.2-6          curl_4.3.2          
##  [91] reprex_2.0.1         bslib_0.3.0          stringi_1.7.4       
##  [94] highr_0.9            ps_1.6.0             Brobdingnag_1.2-6   
##  [97] gdtools_0.2.3        lattice_0.20-44      Matrix_1.3-4        
## [100] nloptr_1.2.2.2       markdown_1.1         shinyjs_2.0.0       
## [103] tensorA_0.36.2       vctrs_0.3.8          pillar_1.6.2        
## [106] lifecycle_1.0.0      jquerylib_0.1.4      bridgesampling_1.1-2
## [109] data.table_1.14.0    httpuv_1.6.3         R6_2.5.1            
## [112] promises_1.2.0.1     renv_0.14.0          gridExtra_2.3       
## [115] codetools_0.2-18     boot_1.3-28          colourpicker_1.1.0  
## [118] MASS_7.3-54          gtools_3.9.2         assertthat_0.2.1    
## [121] rprojroot_2.0.2      withr_2.4.2          shinystan_2.5.0     
## [124] mgcv_1.8-36          parallel_4.0.4       hms_1.1.0           
## [127] grid_4.0.4           coda_0.19-4          minqa_1.2.4         
## [130] cmdstanr_0.4.0       rmarkdown_2.10       shiny_1.6.0         
## [133] lubridate_1.7.10     base64enc_0.1-3      dygraphs_1.1.1.6